今天介紹雷達資料,其實精準地來說,是要看該雷達的功用,才能知道可以得到的資料。
譬如氣象雷達的功用,是透過發射電磁波,當電磁波遇到空氣中有足夠大的雨滴或冰晶時,就能夠測得相關資訊。
而最常用的相關資訊,就是雷達回波。故先介紹讀取雷達回波。
先至氣象局開放資料下載雷達回波資料,如下圖
因為之前已經介紹過如何讀取xml,故下載xml檔案格式,下載之後應該會得到O-A0059-001.xml檔案
使用網頁瀏覽器開啟檔案,如下圖
可以觀察到所需的標記不多,所以就直接iter想要的標記。
處理的方式如下
import xml.etree.ElementTree as ET
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.colors as mcolors
import cartopy.crs as crs
import cartopy.io.shapereader as shpreader
from cartopy.feature import ShapelyFeature
trees = ET.parse('O-A0059-001.xml')
roots = trees.getroot()
pN = "{urn:cwb:gov:tw:cwbcommon:0.1}parameterName"
rN = "{urn:cwb:gov:tw:cwbcommon:0.1}radarName"
pV = "{urn:cwb:gov:tw:cwbcommon:0.1}parameterValue"
parval = []
for hdvalue in roots.iter(pV):
parval.append(hdvalue.text)
dictparm = dict()
cnt = 0
for hd in roots.iter(pN):
if hd.text == '整合雷達名':
for hdvalue in roots.iter(rN):
dictparm.update({hd.text:hdvalue.text})
else:
dictparm.update({hd.text:parval[cnt]})
cnt +=1
print(dictparm)
ac = []
for qq in roots.iter("{urn:cwb:gov:tw:cwbcommon:0.1}content"):
ac = qq.text
reflect = ac.split(",")
reflect = np.array([float(x) for x in reflect]).reshape(881,921)
fig = plt.figure(figsize=(12,8))
axs = plt.axes(projection=crs.PlateCarree())
tw_shp = ShapelyFeature(shpreader.Reader("C:\workplace\ithome\/tw_shp\COUNTY_MOI_1090820.shp").geometries(), crs.PlateCarree(),facecolor="none",edgecolor='k')
lon,lat = np.meshgrid(np.linspace(115.0, 115.0+0.0125*920,921), np.linspace(18.0, 18.0+0.0125*880, 881))
cwb_data = ["#F2F2F2","white","#07FDFD","#0695FD", "#0203F9","#00FF00","#00C800","#019500","#FEFD02","#FEC801","#FD7A00","#FB0100","#C70100","#950100","#FA03FA","#9800F6"]
cmaps = mcolors.ListedColormap(cwb_data,'radar')
colorlevel = [-1000, -900, 0 ,5,10,15,20,25,30,35,40,45,50,55,60,65,70]
norms = mcolors.BoundaryNorm(colorlevel, cmaps.N)
axs.gridlines(draw_labels=True)
pcm = axs.pcolormesh(lon,lat,reflect,cmap = cmaps, norm=norms,transform=crs.PlateCarree())
axs.add_feature(tw_shp)
axs.set_extent([lon[0,0],lon[-1,-1],lat[0,0],lat[-1,-1]])
plt.title(dictparm["時間"])
#故意做一個colorbar
cwb_data = ["#07FDFD","#0695FD", "#0203F9","#00FF00","#00C800","#019500","#FEFD02","#FEC801","#FD7A00","#FB0100","#C70100","#950100","#FA03FA","#9800F6"]
colorlevel = [ 0 ,5,10,15,20,25,30,35,40,45,50,55,60,65,70]
cmaps = mcolors.ListedColormap(cwb_data,'radar')
norms = mcolors.BoundaryNorm(colorlevel, cmaps.N)
plt.colorbar(mpl.cm.ScalarMappable(norm=norms, cmap=cmaps),pad=0.1)
這樣就可以讀取雷達回波跟視覺化囉
視覺化的圖如下
驗證一下,取氣象局相同時間的雷達回波圖
看起來應該是對的喔~